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Abstract 

It is a classical theorem of Liouville that Hamiltonian systems preserve volume in phase space. Any sym- 
plectic Runge-Kutta method will respect this property for such systems, but it has been shown that no 
B-Series method can be volume preserving for all volume preserving vector fields (BIT 47 (2007) 351-378 
& IMA J. Numer. Anal. 27 (2007) 381-405). In this paper we show that despite this result, symplectic 
Runge-Kutta methods can be volume preserving for a much larger class of vector fields than Hamiltonian 
systems, and discuss how some Runge-Kutta methods can preserve a modified measure exactly. 
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1. Introduction 

The construction of numerical schemes for solving ordinary differential equations (ODEs) such that 
some qualitative geometrical property of the analytical solution is preserved exactly by the numerical solu¬ 
tion is an area of great interest and active research today as part of the field of Geometric Integration. The 
most developed topic in this context is that of integrating Hamiltonian systems while preserving the sym- 
plecticity of the flow, and it was found that a class of Runge-Kutta (RK) methods, now called symplectic 
Runge-Kutta (SRK) methods, provides a convenient way to achieve this [ilj, §VI.4]. 

It is a classical theorem due to Liouville that Hamiltonian systems are also volume preserving: for all 
bounded open sets Q. of phase space, the flow map tp t satisfies vol((£,(Q)) = vol(Q) for all t. Equivalently, the 
Jacobian determinant, det(^'(.*)), is 1 for all x and t [jl], VI.9]. Any symplectic mapping of phase space has 
this property, and therefore SRK methods are volume preserving for Hamiltonian systems. Beyond Hamil¬ 
tonian systems, an ODE x = f(x) is volume preserving if and only if / is divergence free (sometimes called 
source free). General volume preservation like this can be found in applications involving incompress¬ 
ible fluid flows and vorticities, ergodic theory and statistical mechanics, and problems in electromagnetism 

Ml- 

One can ask if any SRK methods are volume preserving for all divergence free vector fields /, and 
it has been known for 20 years that the answer is no. Kang and Zai-Jiu showed that no RK method can 
be volume preserving even for the class of linear divergence free vector fields 0 . It was later shown by 
Iserles, Quispel and Tse and independently by Chartier and Murua that no B-Series method can be volume 
preserving for all divergence free vector fields 00 . However, Hairer, Lubich and Wanner have considered 
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separable divergence free vector fields of the form 


fix, y ) = (u(y), v(x)) T , (HLW) 

for functions u : R" — > R m , v : R m — > R" [Q], Thm. 9.4]. There the authors prove that any SRK method with 
at most two stages (and by the product rule of differentiation any composition of such methods) is volume 
preserving for these systems, giving a hint at the fact that SRK methods can be volume preserving for a 
much larger class of vector fields than just Hamiltonian systems. 

As we will show in the introduction, vector fields / in that class must satisfy the determinant condition 

det |/ + -/'(a')| = det |/ - —f(x) j for all h > 0, x e R”, (det) 

where / denotes the nXn identity matrix. In order to substantiate this claim and in anticipation of some of 
the results to be discussed later, we consider the following three Runge-Kutta methods x d>h(x) which 
have been shown to preserve certain measures n(x)dx for quadratic Hamiltonian vector fields |6j]: 

1. The implicit midpoint rule 


<t>h(x) - x _ j ^phix) + v j 


with fi(x) - 1, 


2. the trapezoidal rule 

—- = \{fix) + f (<f>h(x)j), with n(x) = det(l - §/'(*)), (1.1) 

3. and Kahan’s method (restricted to quadratic vector fields) 

—X _ + - jf(x) - \fi<phix)), with nix) = det(l - %f{xj) (1.2) 

These quadratic Hamiltonian vector fields satisfy the determinant condition (Idetl ) and we will establish in 
section |4] that this condition is essential for these measure preservation properties. Indeed, using the chain 
rule, we compute the Jacobian matrix of the midpoint rule to 

which in turn gives the condition for volume preservation 

det (/ + jfi(x + <f> h (x))/ 2)) 

det(0,(x)) =- - -= 1. 

det(7 - |/'((x + Mx))/2)) 

Note that in agreement with |0], it is clear that for the implicit midpoint rule we cannot consider a class 
of vector fields any larger than this and realistically expect volume preservation. Hence we restrict our 
discussion to vector fields satisfying this determinant condition (Idetl) . These functions, as we show later, are 
divergence free and include Hamiltonian systems and HLW separable systems described above. 
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The contributions of this paper are to highlight the relevance of the determinant condition (Idetl ) for 
volume preservation by Runge-Kutta methods, and to introduce and prove results regarding volume preser¬ 
vation for some classes of vector fields lying between Hamiltonian vector fields and those satisfying the 
determinant condition (Idetl ). Not only does this further the understanding of Runge-Kutta methods and vol¬ 
ume preservation of numerical methods in general, but it gives examples of where in applications one could 
in principle use Runge-Kutta methods and preserve volume for a non-Hamiltonian system. Furthermore, 
we discuss how Runge-Kutta methods can also preserve a modified measure exactly. The importance of 
such methods is that the dynamics of the numerical solution lie in the class of measure preserving systems, 
giving a qualitative advantage over methods lacking this property OD. It should be noted that there are gen¬ 
eral approaches to constructing volume preserving splitting methods for a general divergence free vector 
held 01,0,01, but Runge-Kutta methods offer practical and theoretical simplicity and familiarity. 

2. Properties of Runge-Kutta methods 

This section is fairly technical, but it provides us with the necessary tools for the discussion in sections 
3 and 4. We use the following notation to describe a Runge-Kutta method for the autonomous system 
x = f(x). We assume / is continuously differentiable throughout the paper. For each step-size h, a .v-stagc 
Runge-Kutta method provides a map <f>/, : M" —> M”, defined by 

<f>h(x) = x + h^ 

(=t 


where the stages kj satisfy 

kj = x + h ^ cijjf(kj), for i = 1 , ..., s. 
j= i 

As usual, we consolidate the bj s and cij/s into the Butcher tableau consisting of the vector b and the matrix 
A. We make use of the Kronecker product throughout, which for A e M" x '' and B e M' nxm is defined to be 


A®B = 


flu B 


Bni B 


ci in B 


(Inn By 


Lemma 2.1. The Jacobian matrix of a RK method can be written as 

<p' h (x) = 1 + h{b T <8> I)F(I S 0 7- h(A <g> 7)F) _1 (1 0 7), 

with determinant 


det (<p' h (x)) = 


det(7, 0 / - h{{A - lb T ) 0 1)F) 


det(7, 0 7- h(A 0 1)F) 

where F = diag(/ r (ki),... ,f'(k s )), 1 is an s X 1 vector of l’s and I s is the s X s identity matrix. 
Proof Computing directly, we find 

f’ h {x) = 7 + h'Y j b i f\kj)k' j (x) = I + h(b T ® I)F(k[(x),k' s (x)) T . 

i= 1 


( 2 . 1 ) 

( 2 . 2 ) 


(2.3) 
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By definition of the stages ku the derivatives k'(x) satisfy 


7 - ha n f'(ki) 

-hauf'(k 2 ) ■ 

- ha\sf(k s ) N 

(k\(x)\ 


r 

-ha 2 \f(ki ) 

I - ha 2 2 f'(k 2 ) ■ 

- ha 2s f'(k s ) 

k' 2 (x) 

- 

i 

-ha s \f(k\) 

~ha s2 f(k 2 ) • 

■ 7 - hassf'(ks ), 

X s (x), 


J , 


Written more compactly using Kronecker products, this is 

(I s 0 I - h(A ® I)F)(k\ (x),k' s (x)) T = 1 0 7. (2.4) 

The form of the Jacobian matrix can now be found by substituting (12.41) into (12.31 ). 

For the determinant, use the block determinant identity 

det([/)det(X - WU~ l V) = det|^ ^ = det(X)det(U - VX~ l W) (2.5) 

on the expression (12.11) with U = 7 v 0 1 - h(A 0 1)F, V = (1 0 7), W = -h(b T 0 1)F and X = 7. □ 

We wish to understand for which vector fields / and which Runge-Kutta methods defined by A and b, 
the determinant (12.21 ) is unity. As one might expect, this turns out to be simpler for symplectic Runge-Kutta 
methods. Now, for the puipose of exposition, we restrict to methods described in the following definition 
and instruct the reader in how certain results can be proven for general SRK methods at the end of the 
section. 

Definition 2.2. A SRK method is said to be a special symplectic Runge-Kutta method (SSRK) if bj + 0 
for all j, so that the Butcher tableau may be written A = ^(O + 11 T )7?, where B = diag(/;) and Q is a 
skew-symmetric matrix. 

This definition is reasonable because if bj + 0 for all j, then the matrix M = BA + A ' B - bb is zero 
(which implies the method is symplectic) if and only if Q is skew-symmetric. The expression JfiQ + 11 T )B 
therefore constitutes a normal form for most SRK methods of interest [1]. 

Lemma 2.3. An s-stage SSRK method is volume preserving for x = fix) if and only if 

det (I s ®7- h(A ® I)F) = det(7, 0 7 + h(A T ® I)F), (2.6) 

where F = diag(/'(ki),.. .,f(k s )). 

Proof The equation M = 0 can be written -A T = B(A - tb l )B 1 . Hence 

det(7 i ®7 + h(A r ® I)F) = det(7. s 0 7 - h(B(A - tb T )B~ x ® I)F) 

= det(7, ® 7 - h(B 0 7)(A - 1 b T ) ® I)F(B ® 7) _1 ) 

= det(7. s 0 7 - h((A - 1 b r ) 0 1)F) 

The result now follows from Lemma [2Tl □ 
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When s = 1, the only SSRK method is the implicit midpoint rule. In this case. Lemma 1231 gives the 
determinant condition (Idetl) from the introduction. 

When s = 2, we have a three-parameter family of SSRK methods, which reduces to two-parameter if 
we impose the consistency condition b\ + b 2 = 1. Now Lemma [231 gives the condition 

tl-hanfih) -hanfih) \ , (I + hanfih) ha 2 ifih) \ 

\ - ha 21 f\k l ) I - ha 22 f\k 2 )j Qet \ hanfih) 7 + ha 22 f'(k 2 )j ' 

Applying the block determinant identity (12.51) , this boils down to 


det(7 - hanfih) - ha 22 f'(k 2 ) + h 2 det(A)/'(ki)/'(k 2 )) 

= det(7 + hanfih) + ha 22 f(k 2 ) + /; 2 det(A)/'(ki)/'(k 2 )). (2.8) 


We were able here to simplify the identity (12.51) because the top-left block ( I-hanfih )) a nd the bottom-left 
block {-ha 2 \f'{k\)) commute. This cannot be done for ,v > 3. 

These next three lemmata give some basic operations that can be performed on the vector field which 
send volume preserving ODEs to volume preserving ODEs, and effect a simple change in the Jacobian 
determinant of some RK methods for general vector fields. 

Lemma 2.4. Let f : M" —* M" and define a linear change of variables fix) = Pf(P~ l x)for some invertible 
matrix P. Then the RK map fi, for solving x = fix) satisfies 

h(x) = Pfh{P~ X x), ~f h ix) = Pf h iP- l x)p-\ (2.9) 


Lemma 2.5. Let u : M m -4 W n , v : R n+m -» R" and define f : W l,+n -4 R m+n by 


fix,y) 


I u{x) \ 
\v(x,y)j' 


( 2 . 10 ) 


Now let <ph : M" +m —» M' !+m be a one-stage RK map for solving (i,y) T = f(x,y), </p, : M m —> that for 

solving x = u(x), andxn '■ K" +,,i —> R” that for solving y = v(x, y) where x is treated as a parameter. Then 


4>h(x,y) 


'J'hix) \ 
■h(k](x),y)J ’ 


( 2 . 11 ) 


and consequently 

det if h ix,y)) = detif h ix))detid y Xhihix),y)), 

where k\(x) = x + ha\ \ u(k\(x)) is the internal stage of the RK method fifx) and d y denotes the derivative 
with respect to the y coordinate. 


Proof. The full method is 


4>h(x,y) 



+ hb 1 


/ u(k\ (xj) 
\vihix),hihix),y)) 


with the internal stages 


/ hix) \ (x\ , ( uihix)) \ 

\h(h(x),y)J an \vihix),hihix),y))j' 


( 2 . 12 ) 
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The methods applied to each component of the M." +m dimensional system are given by 


/ fifix) 


\Yh(x,y) \y 


+ lib i 


ufkfix)) 
v(x,lfix,y)))' 


(2.13) 


Comparing (12.121) with (12.131) yields the result (12.111) . To prove the last part, note that the Jacobian matrix 
has block structure 

d > ' (x y) = l W ° ) 

h ’ \d x (Xh(ki(x),y)) d y (xh(ki(x),y))J 

and so the determinant det {(p' h {x,y)) is the product of the determinants of the diagonal blocks. □ 

For some simple vector fields, this can be generalized to certain 5-stage methods. Note that the notation 
for Xh is different to that for Lemma [231 


Lemma 2.6. Let u : M"' —> M"', v : W’ —> M' ! , w : M m 


l n and define f ■ M m+ " -» K ' n+n by 


f(x,y) = 


u{x) 
w(x) + v(y) r 


(2.14) 


Now let cp h (x, y) be the RK map for solving (i\y) T = f{x,y), that for solving x = u(x), and Xh{c,y) 

that for solving y = c + v(y). Define c t = a ij■ If the Butcher tableau is such that 


8j(i, k ) = — -— 1 < i, j, k < s, 

Ci - c k 


is finite and independent of distinct i and kfor each j then there exist functions <://,, <?/,, c/, 
that 

fih(x) \ 

Ch(dh(x),y + he h (x)) + he fix) j 


<Pfix,y) = 


for all y. 


Consequently, 


det (<f>' h (x,y)) = det(i/r'(x)) det (d y Xh(dfix),y + he fix))). 


(2.15) 
l n such 

(2.16) 

(2.17) 


Proof Write <pfix,y) = (fih(x), crfix.y)). Note that crfix,y) + Xhfifix),}’), but they are related as follows. 

S S 

crfix,y) - y + h ^ b-Mk,) + h ^ /j,v(/,(vv(/ci), .. ,,w(k s ),y)), (2.18) 

(=t i=t 

S S 

Xh(c,y) = y + h E b ‘ c+h Tj biv(li(c,...,c,y)). 


i= 1 


i= 1 


with stage values 


S S 

Iff l, • • •, fs,y) = y + h ^ aijfj + h ^ ciijvdfifi ,.. . ,f s ,y)). 
j= i 7=1 

Now let d be an arbitrary number. Then we have for each i, 

S S 

h(w(k ]),..., w(k s ), y) = y + lie ,• + h ^ aijd + li ^ a tj v(l fiw{ki ),..., w(k s ), y)), 

7=1 7=1 
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(2.19) 


where e t = Y=\ aij(w(kj) - d). Hence 

lj(w(k\),..., w(k s ),y ) = Ifd, ...,d,y + liej). 


We want to choose d such that <?, = e k Vi, k. Equivalently, 

S S 

Y, a.ij(w(kj) - cl) - Yj a k j(w(k j) - d) for all i + k. 

7=1 7=1 

Solving for d we find 

w(kj) ( — -—) for all i + k. 

\ Ci~C k j 

This will only give us a unique finite value of d no matter what values w(kj) take if the value of 6j(i, k) 
is finite and independent of distinct i and k for every j, which is given by assumption. Hence we can set 
dh(x) = d, e/fx) - e\ and by (12.191 ). we write (12.181) as 



c Th(x,y ) = y + h Y bjw(ki) + li ^ biv(li(d h ,..., d h ,y + he h )) 

i= 1 *=1 

s s s 

= (y + he h ) + h ^ bjd h + li ^ biv(li(d h ,d h ,y + he h )) + h ^ bi(w(ki) - d h ) - 


eh 


i= 1 i=l 

-Xh(d h (x),y + he h (x)) + hc h {x). 


V i=l 


where ci,(x) = (H - =1 b^wiki) - dh) - e/,j. The factorisation of the determinant is evident from the block 
structure of the Jacobian matrix 


<P' h (-x,y) 


Y'h (x ) 0 

\ * d y (xh(dh(x),y + he h (x)) + hc h (x)) 


□ 


Remark 2.7. Let us shed some light on the meaning of (12.151) being finite and independent of distinct i and k 
for each j. The finiteness implies that the method has c -, + cy for all i + k, which is known as nonconfluency 
|[0]. For two-stage SSRK methods, condition (12.151) is satisfied if the method is consistent and f > + 0. There 
is a one-parameter family of self-adjoint three-stage SSRK methods of order four that satisfy the condition. 
The three-stage Gauss-Legendre method, however, does not belong to this class. 


Definition 2.8. A vector field / : M” +m —> W ,+m possesses a linear foliation if there exists a linear change 
of variables as in Lemma 1241 such that / is as in (12.101) from Lemma [231 for some functions u and v. Such 
vector fields are called linearly foliate. See [ 81 ] for general Lie group foliations in the context of Geometric 
Integration. 


Remark 2.9. For general SRK methods, the condition in Lemma [231 along with the condition with A re¬ 
placed by A - 1 b T is sufficient for volume preservation. This result can be obtained along the lines of 
|yj, Thm. 9.4] regarding separable systems (HLW), as follows. Consider the foliation x = fix), y = 
—f'(x) T y, which is Hamiltonian with respect to H(x,y) = y T f(x). Then, using the notation of Lemma 

12.61 the Jacobian matrix of the Runge-Kutta map has block structure ( ^ ^ ^ . |. As in [ 2 , 

\ * dyerh(x, y) J 
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Thm. 9.4], since the vector field is Hamiltonian, a SRK method will produce a symplectic map, which im¬ 
plies det(^(x))det(d y 07 ,(x,;y)) = 1. Hence to show that dct(</>' ; (x)) = I it suffices to show that det(<^(x)) = 
det (d y CTh{x,y)). Computing these two sides as in Lemma 12.11 using the block determinant relation and 
equating numerators and denominators, gives the 2 conditions mentioned above. 

3. Classification of volume preserving vector fields 

Definition 3.1. Define the following classes of vector fields on Euclidean space recursively using vector 
fields f{x,y) = (w(x), v(x, y)) T possessing linear foliations as in Definition 12.81 

P7 = {/ such that there exists P such that for all x, Pf'(x)P 1 = -/'(x) T } , 

S = {/ such that there exists P such that for all x, Pf'(x)P 1 = —f(x)\, 

T ica) = \f(x, y) = (u(x), v(x,y)) T where u e 77 U 'f <03> and there exists P such that for all x,y 
Pd y v(x,y)P ,_1 = -d v v(x,y) T }, 

T (T) = lf(x,y) = (u(x), v(x, y)) T where u e S U ‘H U T (2) and there exists P such that for all x,y 
either Pd y v(x,y)P~ 1 = -d y v(x,y) T or Pd y v(x,y)P~ l = -d v v(x,y)}, 

fD = | vector fields satisfying det(7 + -f{x)) = det(7 - — f'(x)) for all h > 0 and all xj . 

Lemma 3.2. The set '77 contains all vector fields of the form /(x) = 7 _1 V77(x) where J is constant and 
skew-symmetric. All SRK methods are volume preserving for vector fields in '77. 

Proof. For the first part, note that if /(x) = 7 _l V77(x), then Jf'(x)P 1 = V 2 H(x)J 1 = -/'(x) T . For the 
second part, let A e R vXi and P be such that for all x, Pf\x)P~ x = -/'(x) T . Then using the notation of 
Lemma 1231 

det(7, ® 7 - h(A ® I)F) = det(7, ® 7 - h(I s ® P)(A ® 7)(7, ® P _1 )(7, ® P)F(I S ® P _1 )) (3.1) 


= det(7, ® 7 + h(A ® 7)P T ) (3.2) 

= det(7 v ® 7 + hF(A T ® 7) (transpose) (3.3) 

= det(7 v ® 7 + h(A r ® 7)P) (Sylvester’s law). (3.4) 

By Lemma [23] and Remark [23] all SRK methods are volume preserving. □ 


Lemma 3.3. The set S contains all separable HLW systems i.e. fix, y) = (u(y), v(x)) T . All SRK methods 
with at most 2 stages, and compositions thereof are volume preserving for vector fields in S. 

Proof. For the first part, note that if f{x,y) = (u(y), v(x)) T . then Df'(x,y)D 1 = - f'(x,y) where D = 
diag(7 m , -I n )- For the second part, let A e M 2x2 and P be such that for all x P/'(x)P _1 = - fix). Then for 
the two stages k\, ki of the SRK method, 

det(7 - ha\\ f(k\) - ha^f'fa) + h 2 detiA)f(k 1 )f(k 2 )) 

= det(7 - ha\\Pf {k\)P~ x - ha 22 P.f (k 2 )P~ l + h 2 det(A)P/ , (ki)p- 1 P/ , (k 2 )P _1 ) 

= det(7 + hanfih) + ha 22 f(k 2 ) + h 2 det( A)f\k ] )f\k 2 )). 



By (12.81 ) and Remark [2T9l all 2-stage SRK methods are volume preserving. To complete the proof, note that 
a 1-stage SRK method is equivalent to a 2-stage SRK method with two equal stages, and compositions of 
volume preserving maps arc also volume preserving. □ 

Lemma 3.4. The inclusions 77 c (T' (oo) c T c 2) and S c T (2) c 1) hold. 


Proof. PI c p' {ay) c T {2) and S c T {2) are clear by considering trivial foliations in which n + m = m. We 
will show that S c 2), PI c 2) and that 1) is closed under the employed recursive process leading to linearly 
foliate systems. 

For / € S , det(7 + 'ff’(x)) = det(7 + \Pf'(x)P~ l ) = det(7 - §/'(x)). 

For / 6 *77, det(7 + §/'(x)) = det(7 + \Pf(x)P~ x ) = det(7 -~§/'(x) T ) = det(7 - §/'(*)). 

Let / £ 2) and define /(x) = Pf(P~ l x ) for an invertible matrix P. Then det(7 + j f'(x)) = det(7 + 
| Pf (P~ l x)P~ x ) = det(7 + ^f(P' 1 v). Doing the same with a - instead of a + shows that / € £). 

Let f{x,y) = (u(x), v(x, y)) T where u e D and y v(x,y) € 2) for all v. Then 


det ( 7 + 


(I + ju'(x) 0 \ 

Ct \ \d x v{x, y) I + \d y v(x, y)j 

det(7 + -^u'(x)) det(7 + ^d y v(x,y)). 


(3.5) 

(3.6) 


Doing the same with a - instead of a + shows that / £ 2). 


□ 


Theorem 3.5. The following are equivalent. 

(i) f € £> 

(ii) det(7 + zf'(x)) = det(7 - zf\x))for all z £ C tmr/ fl// x 

(Hi) The non-zero eigenvalues of f'(x), counting multiplicities, come in positive-negative pairs 
(iv) U'( f'(x) 2k+l ) = 0 for all x and k = 0,1,2,... 

Proof. (0 <=> (ii): Assuming (;), for every x, p(z) = det(7 + zf'(x)) - det(7 - zf'(x)) is a polynomial in z 
that is zero for infinitely many values of z = h/2 e M.+. By the Fundamental Theorem of Algebra, p(z) — 0 
for all z £ C. The converse follows from setting h = 2z £ M+ c C. 

(ii) => (Hi): By triangularisation we can see that for every x, the polynomial q(z) = det(7 - zf'(x)) 
is equal to (1 - zdi) •■•(!- z.A r ) where ..., A r are the non-zero eigenvalues of f'(x). If (i) holds, then 
q(z) - q(~z), and the roots I /A,- of q come in positive-negative pairs. Hence the eigenvalues Aj do too. 
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(Hi) ==> (iv): For all x, tr(f'(x) 2k+1 ) = /l 2/i+l -t-b A 2k+ 1 where A\,..., A r are the non-zero eigenvalues 

of f'(x). Flence if the non-zero eigenvalues come in positive-negative pairs then tr (f'(x) 2k+1 ) = 0 for 
k - 0,1,2,.... 

(iv) ==> (it): Newton’s identity gives 

^ 2k+l 

e 2k+ iUi ,...,A n ) = r—r Y(-l t l e 2M -i(A x , ■ ■ ■ ,A n )Xx(f'(xj), (3.7) 

2k + 1 

i= l 

where efA],... ,A n ) is the elementary symmetric polynomial in A\, ... ,A n and, incidentally, the coefficient 
of z) in q(z) = det(7 - zf'(x)). Since for any k and i, either 2k + 1 - i is odd or i is odd, we can use an 
induction argument to show that all the coefficients of z 2k+i in q(z) are zero. Flence q(z) = q(~z)- □ 

Corollary 3.6. All elements of A) are divergence free. Restricted to 2 dimensional vector fields, D, 77 and 
S are all equal to divergence free vector fields. 

Theorem 3.7. The set 'F (oo) contains all 

(i) Affine vector fields f(x) = Lx + d such that det(7 + |L) = det(7 - l fL)for all h > 0 

(ii) Vector fields such that f'(x) = JS (x) where J is skew-symmetric and S (x) is symmetric 


Proof, (i) Let L satisfy the determinant condition (Idetl) . By the Jordan normal form, and the fact that the 
eigenvalues must come in positive-negative pairs by Theorem 13.51 we can find an invertible matrix P such 
that 

PLP~ l = diagfiljf + N u -A^l + N- X ,A 2 1 + N 2 ,-A 2 I + AL 2 ,..., A,I + N r ,-A r I + AL r , N 0 ), 


where the TV* are matrices that are zero everywhere except for possible l’s on the first subdiagonal (Nk)i+ij- 
Hence / is a tower of linear foliations of affine functions with Jacobian matrices either No or diag(T7 + 
N i, -/17 + N-\). If f(x) = Nqx + cl then this is clearly a tower of foliations of zero systems i.e. u(x) = 0, 
v(x,y) = x. Now consider the case f(x) = diag(d7 + N\, -AI + N-\)x + d. There is a simple permutation of 
variables so that the Jacobian matrix becomes 


'A 

0 

0 

0 

0 

0 

0 ' 

0 

-A 

0 

0 

0 

0 

0 

★ 

0 

A 

0 

0 

0 

0 

0 

★ 

0 

-A 

0 

0 

0 

0 

0 


★ 

0 

A 

0 

lo 

0 


0 

★ 

0 

-A, 


where the *’s are possible l’s (0 otherwise). Hence / is a tower of linear foliations of harmonic oscillators, 
u(x ], x 2 ) = (Axi , -Ax 2 ), v(x\,x 2 ,yi,y 2 ) = (*x\ + Ay 1 , *x 2 - /ly 2 ). 

(ii) By a linear orthogonal change of variables, we can assume J = diag(0, 77 1 ), where K is skew- 
symmetric. In this case there is symmetric T(x) and V(x) such that 


fix) = 


0 

K - 1 


I T(x) 
\U(x) T 


U(x)\ 

V(x)j 


0 0 \ 
K~ l U(x) T K~ l V(x)J' 


(3.8) 


This shows that / possesses a linear foliation with a zero system u e *77 and a system v with d y v(x, y) = 
K~ l V(x) so that y i-> v(x,y) € Li with the same P = K for all x,y. □ 
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Theorem 3.8. Consider an s-stage SRK method that is volume preserving for the vector field u : K"' —» W n , 
and let v : W n+n —> W' 1+n be such that there exists an invertible matrix P such that for all x,y, 

Pdyv(x,y)P~ ] - -d y v(x,y) T . 

Then the SRK method is volume preserving for the vector field 

fix, y) = (; u(x ), v(x, y )) T . (3.9) 


Proof Let A e M sx,v and take P from the assumption. By Lemma [231 and Remark 12.91 a SRK method is 
volume preserving if 

det(4 0 1 - h(A 0 1)F) = det(7 s 0 / + h{A r 0 1)F), (1231) 

where F = diag(/'(ki),..., f'(k f ). For (13.91) . the Jacobian matrix becomes 

, lu\x) 0 \ 

X,y \ * d y v(x,y)j 

and using a similarity transformation, we can bring det(7 0 1 - h(A 0 I)F) to the form 


det 



-flljKj 


0 

0 ' 

-a s \ u\ 

■ 7 - a ss u' s 


0 

0 

★ 

★ 

7 

-a n v\ ■ 

• a\ s v' s 

v ★ 

★ 

7 

- a s \v\ ■ 

ClssV s' 


(3.10) 


where it'-, v' are shorthand for d x u(ki) and d Y v{kj), respectively. Thus, the condition (12.61) factorises to 


det(7, 0 / - h(A 0 1)U ) det(/ s 0 / - h(A 0 1)V) = det(7, 0 7 + h{A r 0 7)77) det(7, 07 + h(A T 0 1)V), 


with U = diag(i<',..., u' s ), V = diag(v'j,..., v'). We compute 

det(7, 0 7 - h(A 0 1)V) = det(7, 0 7 - h(I s 0 P)(A 0 7)F(7, 0 P) _1 ) 
= det(7, 0 7 - h(A 0 7)(7, 0 P)V(I S 0F 1 )) 
= det(7, 0 7 + h(A 0 1)V T ) 

= det(7, 0 7 + hV(A T 0 7)) 

= det(7, 0 7 + h(A r 0 7)V). 


The last line comes from Sylvester’s determinant identity. The proof is completed noticing that det(7 s 0 7- 
h(A Cl)U) = dct(/ v 07 + h(A 1 0 7)77) is satisfied by the assumption that the method is volume preserving 
for u. □ 

Corollary 3.9. All SRK methods are volume preserving for vector fields in f(°°\ 

Proof SRK methods are volume preserving for vector fields in FI by Lemma [L2l and volume preservation 
for the recursive constructions of F (ca> is assured by Theorem 13.81 □ 
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Theorem 3.10. Consider a SRK method with at most two stages (or a composition of such methods) that is 
volume preserving for the vector field u : R' n —* R m , and let v : M m+,! —> M.' n+n be such that there exists an 
invertible matrix P such that for all x,y, 

Pdyv(x,y)P~ l = -d y v(x,y). 

Then the SRK method is volume presen’ing for the vector field 

f(x, y) = (u(x),v(x,y)) T . 

Proof. Let A e M 2x2 and take P from assumption. As in Theorem 13. 8 1 the Jacobian matrix is block triangu¬ 
lar 

f'i t = ( “’W 0 \ 

1 X ’ y) \d x v(x,y) d y v(x,y)]- 

For 2-stage methods, the condition for volume preservation from equation (12.81) is 


det(7 - haufXh) - ha 22 f(k 2 ) + h 2 det(A)/'(ki)/'(k 2 )) 

= det(7 + hanffki) + ha 2 2 f(k 2 ) + h 2 dct(A)f(k x )fXk 2 )). 

Now, because of the block-triangular structure of /'(k,) and 


f'(ki)f'(k 2 ) = 


u'(ki)u'(k 2 ) 


0 

d y v(k\)d y v(k 2 ) 


, f(k\) + f'(k 2 ) = 


u'(k\) + u'(k 2 ) 


0 

dy\’(k\) + dyv(k 2 )] ’ 


where we have used the convention that u(kj) has used the x component of k\. The condition (12.81 ) then 
factorises into 


det(7 - h(a n f(ki) + a 22 f'(k 2 )) + h 2 det(A)f(k x )f(k 2 )) 

= det(7 - h(a\\u'(k\) + a 22 u(k 2 )) + h 2 det(A)u'(k x )u'(k 2 )) 

■ det(7 - h(ai\d y v(ki) + a 22 d y v(k 2 )) + lr det(A)d y v(ki)d y v(k 2 )). 

A similarity transformation with P leads to 

det(7 - h{a\\d y v(k\) + a 22 d y v(k 2 )) + lr det(A)i9 ;y v(ki)(9 v v(k2)) 

= det(7 - h(a n Pd y v(kx)p- x + a 22 Pd v v(k 2 )P~') + h 2 det(A)PdyV(k x )P~ l Pd y v(k 2 )P~ l ) 

- (3.11) 

= det(7 + h(a xx dyv(k\) + a 22 d y v(k 2 )) + lr det(A)d y v(k])d y v(k 2 )). 

Condition (12.81 ) for / is now satisfied by considering (12.81 ) for the vector field u (which holds because we 
assume the SRK method is volume preserving) and (13.1 II) . This proves the result for 2-stage SRK methods. 
To complete the proof, note that a 1-stage SRK method is equivalent to a 2-stage SRK method with two 
equal stages, and compositions of volume preserving maps are also volume preserving. □ 

Corollary 3.11. All SRK methods with at most two stages (and compositions thereof) are volume presen’ing 
for vector fields in T {2) . 

Proof. All SRK methods with at most two stages (and compositions thereof) are volume preserving for 
vector fields in dd by Lemma [3721 and S by Lemma 1331 Volume preservation for the recursive construction 
of T {2) is assured by Theorems 13.81 and 13.101 □ 


12 






We already saw in the introduction that the implicit midpoint rule (which is the only 1-stage SRK 
method) is volume preserving for all /efl, and that all such vector fields must lie in D. However, does the 
set T {2) contain all vector fields such that all 2-stage SRK methods are volume preserving? And does the 
set T (ca) contain all vector fields such that all SRK methods are volume preserving? We do not yet know 
the answers to these questions, but the following counterexamples are relevant. 

The first counterexample shows that Corollary 13.1 II is not true for three-stage methods. In the second 
example, we show that D \ r F (2> + 0 and that only the midpoint rule can be volume preserving for all 
methods in T). This counterexample is of the lowest possible dimension (3) but one might argue that 
volume preservation is hindered in this example because the vector field is not completely smooth at x - 0. 
The third example clarifies the matter: we give a way to construct a class of (smooth) vector fields in D for 
which two-stage methods cannot be expected to preserve volume. 

Example 3.12. Hairer, Lubich and Wanner [T , VI.9] used the vector field 


x = sinz, y = cosz, z = siny + cos x, 


to show that the three-stage Gauss-Legendre method is not volume preserving, despite the vector field lying 
in S. What could be interesting is to find some class of functions F <2) such that all three-stage SRK methods 
are volume preserving, but not all four-stage SRK methods. 


Example 3.13. Consider the continuously differentiable vector field 


f(x, y,z) 


(^x 3 - c, -x 2 y, 0) if x > 0 
(^x 3 - c, 0, -x 2 z) if x < 0 



( x 2 

0 

O' 


f x 2 

0 

0 ' 

f'(x,y,z) = 

-2 xy 

-x 2 

0 

if x > 0, 

0 

0 

0 


, O’ 

0 

oj 


-2 xz 

0 



Then / e D, but not all 2-stage SRK methods are volume preserving. The principle here is that if k\ 
and k 2 have x-components with different signs, then f'{k\) and /'(G) will violate the condition in (12.81) . 
For instance, the two-stage Gauss-Legendre method with initial value (1/2,0,0), drift c = 1 and step size 
h =1/2 has stage values with different signs in the x-coordinate and hence, does not preserve volume. 


Example 3.14. The following example illustrates that SRK methods cannot preserve simple systems in T) 
that do not belong to the classes F {oo) or F {2) . Let g e D and let A(x) be skew-symmetric (and invertible) 
and S (y) be symmetric matrices. Then, any vector field with Jacobian matrix 


f'(x,y ) 


(g\x) 0 \ 

\ * Mx)S(y)) 


will satisfy the determinant condition, however, the similarity transform P to yield Pd y f(x,y)P ~ 1 = -d y f(x, y) T 
is now P = A(x )~ 1 and this dependence on x hinders a crucial step in the above proof. For volume preser¬ 
vation of SRK methods, it is thus essential to have a constant transform P for all values of x,y or at least in 
a region of interest for the numerical integration. We give the following concrete example, 



0 X] xi 


O 

O 

CN 

A(x) = 

—xi 0 X1X2 

, S(y) = 

0 y \ 0 


,-xi -X1X2 0 , 


CO 

O 

O 
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which is combined with the harmonic oscillator g(x i, * 2 ) = (xi, -x\) T and could originate from f{x,y) = 
\y\, ^V[) T ) ■ Integrating with step size h = 1/2 from (xo,yo) = (1,1/2,1/3,1/4,1/5) leads 
to a change of volume for the two-stage Gauss-Legendre method. The implicit midpoint rule preserves 
volume as expected. 


4. Measure-preservation of Runge-Kutta methods 


In the introduction, we have pointed out that the trapezoidal method is not necessarily volume preserving 
but instead preserves the measure det (i - §/'(x)) dx for quadratic Hamiltonian vector fields jfij]. Recall that 
a map <f> preserves a measure p{x)dx if 


det (</>'(x)) = 


pjx) 

d(<P(x)) 


This result is generalised in the following lemma. 


Lemma 4.1. The trapezoidal rule (11.11) presents the measure p(x)dx with 

p{x) = det (/ ± | f{xj) 

for vector fields f that satisfy the determinant condition (Idetl) . 


Proof We compute the Jacobian matrix df of the trapezoidal rule, 

<f>h(x) = I + ^/'(x) + ^f'(<f>h(x))<p' h (x). 


and see that 


det .($,(*)) = 


det(7 + f/'(x)) 
det(7 - §/'(<Mx))) 


P(x) 

pi<Ph(x ))' 


□ 


This means that volume is conserved to order 0(h 2 ) globally (by Theorem 13.51) . but more importantly 
that the dynamics of the numerical solution lie in the class of measure preserving systems, giving a qualita¬ 
tive advantage over methods lacking this property (TO. The trapezoidal method is conjugate to the implicit 
midpoint rule 1 1, VI.8], which goes some way towards explaining this behaviour. However, the next method 
we consider has similar measure preserving properties, but doesn’t appear likewise to be “conjugate to vol¬ 
ume preserving”. 

There has been recent interest in the properties of the Kahan method [ji,[9]. For a quadratic vector field 
fix) = Q(x) + L(x) + d where Q is quadratically homogeneous, L is linear and d is constant, the symmetric 
bilinear form q(x,y) is formed by polarisation, 

q(x,y) = ^(<2(x + v) - Q(x) - Qiy )), (4.1) 

and Kahan’s unconventional numerical method is then given by 

—- = qix, (phi x)) + jL (x + </>/,(x)) + d. (4.2) 

In 10], it was shown that Kahan’s method is equivalent to a three-stage Runge-Kutta method restricted to 
quadratic vector fields. We give the following generalisation. 
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Lemma 4.2. Restricted to quadratic vector fields, Kalian’s method is equivalent to the s-stage Runge-Kutta 
method 

S 

4>h(x) = x + h ^ bif{x + Ci(tf>h(x ) - x)), (4.3) 

i=i 

for any b and c satisfying X,' = i bt = 1, X,L i bjCi = X,' = i = 0- This implies that the Butcher tableau 

satisfies A = cb T , but the converse is not true. 

Proof. Let x' = (phlx) and write the vector field as fix) = qlx, x) + Lx + d with the symmetric bilinear form 
q, then, expanding out and setting equal to Kahan’s method (14.21) 

———- = ^ biq{x + Cilx - x), x + Ci{x - xfj + L(x + cfx' - x)) + d 
i= 1 

= q{x', x) + \L{x + x) + d 


yields the above conditions. □ 

In 0 Prop. 5], it was shown that for quadratic Hamiltonian vector fields, Kahan’s method preserves the 
measure with density p(x) - det(7 - ^ f'lx))' 1 . The proof is easily extended to all quadratic vector fields 
satisfying the determinant condition (Idetl ). 


Lemma 4.3. Kahan’s method preserves the measure p{x)dx with 

p(x) = det (/ ± §/'(*)) 1 

for quadratic vector fields f that satisfy the determinant condition (Idetl) . 

Proof. We compute the Jacobian matrix (p' h of Kahan’s method in the form (11.21) . 


flhto = 


I -'jf\x) + 

1 + f/'OfeOO) - \ f ’ 


Since / is quadratic, f is affine and thus 


det {(pjfx)) = 


det(7 + y’Whjx))) 
det(7 - |/'(x)) 


RlsPhlx))' 


(4.4) 


□ 


Due to the similarity of this measure to that preserved by the trapezoidal method, one might at first 
glance suggest that the Kahan method is conjugate to some volume preserving method too, but this does 
not appear to be the case. At least, Kahan’s method is not conjugate by B-series to any symplectic method 
01. It may be interesting to investigate how these measure preserving properties of the trapezoidal rule and 
Kahan’s method can be generalised. 

From Lemmata 12.41 12.51 and 12.61 on linear foliations follow similar measure preservation properties 
generalising the volume preservation properties discussed in the previous section. 

Theorem 4.4. Suppose that a given Runge-Kutta method preserves the measure p on M" when solving the 
ODE x — f{x). Then when solving the ODE x = fix), where fix) = Pf(P~ ] x)for some invertible matrix 
P, the method preserves the measure with density plx) = plP 1 x). 
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Proof. By assumption, det ((p' h (y))p(cf>h(y)) = p(y) for all y e M”. Using the notation and results of Lemma 
El det (<p' h (x))p(<j) h (x)) = det(f’ h (P~ l x))p(f h (P~ x x)) = p(P~'x) = p(x). □ 


Theorem 4.5. Suppose that a given 1-stage Runge-Kutta method preserves the measure pdx on M"' when 
solving the ODE x = u(x), and it preserves the measure v{y)dy on M” when solving the ODE y = v(x, y) for 
all x e M m . Then when solving the ODE (x,y) = (u(x), v(x, y)), the method preserves the product measure 
p(x, y)dxdy = p(x)v(y)dxdy on W l+m . 

Proof. By Lemma [231 <j>h(x,y) - (ifth(x),Xh(ki(x),y)) T , where k\ (x) is the internal stage of the 1-stage 
method. Hence by the definition of p, 

p(<ph(x,y)) = p('l'h(x))v(xh(ki(x),y)). 

By assumption, we have for all jc and y, 

det(f' h (x))p(i// h (x)) = p(x), det(d y Xh(x,y))v(xh(x,y)) = v(y). (4.5) 

Finally, Lemma [231 gives us det{f h {k\{x),y)) = det{f' h {x))det{d y Xh{x,y)). Combining all of these results, 

d&tlf'bix, y))p(<p h (x, y)) = det(f / h (x))p(f, I (x)) det(d y x/i(k, (x), y))v(xh(h (x), y)) 

= P(x)v(y ) 

= p(x,y). 

Hence the measure with density p on is conserved. □ 

From the results of Lemma 12.61 and using its notation, we deduce that a generalization for measure 
preserving RK methods with more stages even for sums f{x,y) = (u(x), w(x) + v(y)) T is not trivial since 
then, 

det(^(x,y)) = det(if/' h (x)) det{d y Xh(d(x),y + he(x))), 
and a product measure p(x, y)dxdy = p{x)viy)dxdy transforms according to 

det (<p' h (x,y))p(<f> h (x,y)) = det (<p' h (x,y))p(ifr h (x)) v(xh(d(x),y + he(x)) + hc(x)) 

= det {f' h {x))p{ip h {x)) ■ det{d y Xh(d(x),y + he(x))} v{xh(d(x),y + he(x)) + hc(x )). 

Assume that fh,Xh preserve the measures with densities p{x) and v(y), respectively, then, if <■/, = <?/, = 0, the 
product measure is preserved. This additional condition holds, e.g., for the trapezoidal rule for which we 
get that dh{x) = (w(k\) + w{kf))l2. Further methods satisfying c>, = eh = 0 can be constructed easihQ but 
they might preserve measures for trivial vector fields only. Kahan’s method derived from Lemma l4~2l does 
not simplify in this way, however, we can give the following result: 

Theorem 4.6. Generalized Kalian’s methods from Lemma \4~2\p reserve the measure p(x, y)dxdy with p(x, y) = 
det(7 + l ff'(x,y))~ ] for linearly foliate vector fields of the form f(x, y) = (u(x), v(y) + w(.r)) T where w is 
arbitrary, and u, v 6 D are quadratic. 


*Let, e.g., aij = 0 and a,- ; - = bj for some ; and all j. 
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Proof. Let z - (x,y), and write fh(z) = {ifh(x), <Th(x,y)) T . We compute the Jacobian determinant of (14.31) 


det ((f>' h (z)) = 


det(7 + h bj( 1 - q)/'(z + c0h(z) - z))) 
det(7 - h biCif'(z + Ci(<p h (z ) - z)))) 


using that /' is block-diagonal, we arrive at 

_ det(7 + h Yjh bi(\ - Ci)u'(x + cffifx) - x ))) det(7 + h Y^ =l bf 1 - c ( )v'(y + a{cr h (x,y) - y))) 
det(7 - /? bjCiu'ix + cffifx) - x ))) det(7 - /? Vff'Cv + Ci{cr h (x, y) - y))) 

and since u', v' are affine, we can simplify using the assumptions on the coefficients from Lemma l4~2l to 

_ det(7 + \u\fifx))) det(7 + \v'{cr h (x,y))) _ det(7 + \ f{(ph(x,y))) _ /u(z) 

det(7 - |n'(x)) det(7 - \v’(y)) ~ det(7 ~'+J\x,y)) ~ piMz))' 

□ 

Remark 4.7. The theorem is not true for more general foliate vector fields within the class p' < °°\ e.g., 
f(x, y) = (u(x), J-' V v 77(v, y)) T where u{x) is a simple harmonic oscillator and with the Hamiltonian H(x, v) = 
iPx ( lx)Py ( ly using the usual notation for the momentum and position coordinates x - ( q x ,Px )> y = ( q y ,p y )■ 
Note that the Hamiltonian is still quadratic in y! 
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